#Motzkin.txt
#This is the Motzkin analogue to Dr. Zeilberger's `Dyck.txt`

print(`To see the procedures in this package, type "Help()"`):
print(`For explanation of each procedure, type "WhatIs(procedure_name)"`):
print(`This package also contains procedures from Dr. Doron Zeilberger's FindRec.txt`):
print(`For a list of the procedures from FindRec.txt`):
print(` type "ezraFindRec();". For specific help type "ezraFindRec(procedure_name);" `):
print(`This package also contains procedures from Dr. Doron Zeilberger's SCHUTZENBERGER.txt`):


Help:=proc(): print(`SetU(m,n)`,`SetD(m,n)`,`SetF(m,n)`): print(`u(m,n)`,`d(m,n)`,`f(m,n)`): 
  print(`Ru(m,n,A,B,F)`, `Rd(m,n,A,B,F)`, `Rf(m,n,A,B,F)`):  
  print(`NRu(m,n,A,B,C,D,E,F)`, `NRd(m,n,A,B,C,D,E,F)`, `NRf(m,n,A,B,C,D,E,F)`):
  print(`NRur(m,n,A,B,C,D,E,F,r)`, `NRdr(m,n,A,B,C,D,E,F,r)`, `NRfr(m,n,A,B,C,D,E,F,r)`):
  print(`SeqABCDE(A,B,C,D,E,K)`, `SeqABCDEr(A,B,C,D,E,r,K)`): 
  print(`Motzkin(n)`, `Mk(n,k)`, `Targem1(L)`, `IsGoodCDE(L,C,D,E)`):
  print(`Targem2(L)`, `IsGoodA(L,A)`, `IsGoodB(L,B)`, `MotzkinABCDE(A,B,C,D,E,n)`, `SeqABCDEbrute(A,B,C,D,E,K)`):
  print(`IsGoodCDEr(L,C,D,E,r)`, `IsGoodAr(L,A,r)`, `IsGoodBr(L,A,r)`, `MotzkinPathsABCDEr(A,B,C,D,E,r,n)`, `SeqABCDErBrute(A,B,C,D,E,r,K)`):
  print(`The following procedures from Dr. Doron Zeilberger's Dyck.txt are also included: IsRange1(f,k,n), IsRange(F,k,n), Aeq(gu,x,P)`):
end:

WhatIs:=proc(procedure):
  if procedure=SetU then
    print(`SetU(m,n): inputs non-negative integers m,n`):
    print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.`):

  elif procedure=SetD then
    print(`SetD(m,n): inputs non-negative integers m,n`):
    print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):

  elif procedure=SetF then
    print(`SetF(m,n): inputs non-negative integers m,n`):
    print(`outputs set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):

  elif procedure=u then
    print(`u(m,n): inputs non-negative integers m,n`):
    print(`outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.`):

  elif procedure=d then
    print(`d(m,n): inputs non-negative integers m,n`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.`):

  elif procedure=f then
    print(`f(m,n): inputs non-negative integers m,n`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step.`):

  elif procedure=Ru then
    print(`Ru(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
    print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
    print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=Rd then
    print(`Rd(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
    print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
    print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=Rf then
    print(`Rf(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f`):
    print(`Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
    print(`s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=NRu then
    print(`NRu(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or d or f`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
    print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
    print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=NRd then
    print(`NRd(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or u or f`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
    print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
    print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=NRf then
    print(`NRf(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or u or d`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
    print(`s.t. no peak has height belongs to A and no valley height belongs to B`):
    print(`no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):
 
  elif procedure=IsRange1 then
    print(`IsRange1(f,k,n): inputs a linear expression f in k and a non-negative integer n`):
    print(`decides whether f(k0)=n for some non-negative integer k0.`):
    print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):

  elif procedure=IsRange then
    print(`IsRange(F,k,n): inputs a set of linear expressions F={f} in k and a non-negative integer n`):
    print(`decides whether f(k0)=n for some non-negative integer k0 and some f in F`):
    print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):

  elif procedure=NRur then
    print(`NRur(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or d or f`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step`):
    print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
    print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=NRdr then
    print(`NRdr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or f`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step`):
    print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
    print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procedure=NRfr then
    print(`NRfr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or d`):
    print(`Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step`):
    print(`s.t. no peak has height is in the range of A and no valley height is in the range of B`):
    print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk.`):

  elif procecedure=SeqABCDE then
    print(`SeqABCDE(A,B,C,D,E,K): inputs finite sets A,B,C,D,E of non-negative integers and positive integer K`):
    print(`Outputs the first K terms of the sequence of Motzkin paths of length n`):
    print(`where no peak height belongs to A and no valley height belongs to B`):
    print(`No upward run belongs to C, No downward run in D, and no flat run belongs to E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk`):

  elif procedure=SeqABCDEr then
    print(`SeqABCDEr(A,B,C,D,E,r,K): inputs sets A,B,C,D,E of linear expressions in the variable r and positive integer K`):
    print(`Outputs the first K terms of the sequence of Motzkin paths of length n`):
    print(`where no peak has height is in the range of A and no valley height is in the range of B`):
    print(`no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E`):
    print(`Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer`):
    print(`Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer`):
    print(`or a chain D F^k that occurs at the end of a walk`):

  elif procedure=Motzkin then
    print(`Motzkin(n): outputs all Motzkin paths of length n.`):

  elif procedure=Targem1 then
    print(`Targem1(L): translates a Motzkin path into a sequence of positive, negative integers, and sets`):
    print(`representing the upwards, downwards, and flat runs respectively,`):
    print(`such that the partial sums of the integers are always non-negative`):

  elif procedure=IsGoodCDE then
    print(`IsGoodCDE(L,C,D,E): Inputs a Motzkin path L, and finite sets C,D,E of non-negative integers.`):
    print(`Decides whether no upwards run in L belongs to C, no downwards run in L belongs to D, and no flat run in L belongs to E.`):

  elif procedure=IsGoodCDEr then
    print(`IsGoodCDEr(L,C,D,E,r): Inputs a Motzkin path L, and sets C,D,E of linear expressions in the variable r.`):
    print(`Decides whether no upwards run in L belongs to the range of C, no downward run in L belongs to the range of D, and no flat run in L belongs to the range of E.`):

  elif procedure=Targem2 then
    print(`Targem2(L): translates a Motzkin path into a sequence of the partial sums when it between a`):
    print(`chain with steps in {0,1} starting with 1 (i.e. weakly increasing chain starting with an upwards step)`):
    print(`and a chain with steps in {0,-1} starting with -1 (i.e. a weakly decreasing chain starting with a downwards step)`):

  elif procedure=IsGoodA then
    print(`IsGoodA(L,A): Inputs a Motzkin path L, and a finite set A of non-negative integers.`):
    print(`Decides whether none of the heights in a peak location of L belongs to A`):
    print(`A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1, s.t. any steps between (if any) are 0.`):

  elif procedure=IsGoodB then
    print(`IsGoodB(L,B): Inputs a Motzkin path L, and a finite set B of non-negative integers.`):
    print(`Decides whether none of the heights in a valley location of L belongs to B.`):
    print(`A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1 s.t. any steps between (if any) are 0.`):

  elif procedure=IsGoodAr then
    print(`IsGoodAr(L,A,r): Inputs a Motzkin path L, and a set A of linear expressions in the variable r.`):
    print(`Decides whether none of the heights in a peak location belongs to the range of A`):
    print(`A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1, s.t. any steps between (if any) are 0.`):

  elif procedure=IsGoodBr then
    print(`IsGoodBr(L,A,r): Inputs a Motzkin path L, and a set B of linear expressions in the variable r.`):
    print(`Decides whether none of the heights in a valley location belongs to the range of B`):
    print(`A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1 s.t. any steps between (if any) are 0.`):

  elif procedure=MotzkinABCDE then
    print(`MotzkinABCDE(A,B,C,D,E,n): Inputs five finite sets A, B, C,D,E of non-negative integers and a positive integer n.`):
    print(`Outputs the set of Motzkin paths counted by SeqABCDE(A,B,C,D,E,K). In other words:`):
    print(`where no peak can be of height that belongs to A and no valley can be of height that belongs to B`):
    print(`No length of upward run belongs to C and no length of downward run belongs to D and no length of flat run belongs to E`):
    print(`A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down), s.t. any steps between (if any) are 0 (i.e. flat).`):

  elif procedure=MotzkinABCDEr then
    print(`MotzkinABCDEr(A,B,C,D,E,r,n): Inputs five sets  A, B, C,D,E of linear expressions in the variable r, and a positive integer n.`):
    print(`Outputs the set of Motzkin paths counted by SeqABCDEr(A,B,C,D,E,r,K). In other words:`):
    print(`where no peak can be of height that belongs to the range of A and no valley can be of height that belongs to  the range of B,`):
    print(`No length of upward run belongs to the range of C and no length of downward run belongs to the range of D and no length of flat run belongs to the range of E.`):
    print(`A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down), s.t. any steps between (if any) are 0 (i.e. flat).`):

  elif procedure=SeqABCDEbrute then
    print(`SeqABCDEbrute(A,B,C,D,E,K): Same output as SeqABCDE(A,B,C,D,E,K) but by complete brute force.`):
    print(`FOR CHECKING PURPOSES ONLY; Don't  make K too big`):

  elif procedure=SeqABCDErBrute then
    print(`SeqABCDErBrute(A,B,C,D,K):  Same output as SeqABCDEr(A,B,C,D,E,r,K) but by complete brute force.`):
    print(`FOR CHECKING PURPOSES ONLY; Don't  make K too big.`): 
 
  elif procedure=Aeq then
    print(`Aeq(L,x,P): The guessed algebraic equation in P(x) satisfied by the generating function of the sequence whose beginning`):
    print(`starting at n=0 is L. Try:`):
    print(`Aeq([seq(binomial(2*i,i)/(i+1),i=1..10)],x,P);`):
    print(`This was taken from Dr. Doron Zeilberger's Dyck.txt`):
  fi:
end:


#SetU(m,n): inputs non-negative integers m,n
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.
SetU:=proc(m,n) local gu,mu,r,mu1 :
  option remember:

  if (m<0 or n<0 or m<n) then RETURN({}): fi:

  if m=0 and n=0 then RETURN({}): fi:

  gu:={}:

  for r from 1 to n do 
    mu:=SetD(m-r,n-r) union SetF(m-r,n-r):
    gu:= gu union {seq([op(mu1),1$r],mu1 in mu)}: 
  od:

  gu:

  end:


#SetD(m,n): inputs non-negative integers m,n
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.
SetD:=proc(m,n) local gu,mu,r,mu1 :
 option remember:

  if (m<0 or n<0 or m<n) then RETURN({}): fi:

  if m=0 then 
    if n=0 then RETURN({[]}):
    else RETURN({}): fi: fi:

  gu:={}:

  for r from 1 to m do 
    mu:=SetU(m-r,n+r) union SetF(m-r,n+r):
    gu:= gu union {seq([op(mu1),-1$r],mu1 in mu)}: 
  od:

  gu:

  end:


#SetF(m,n): inputs non-negative integers m,n
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step.
SetF:=proc(m,n) local gu,mu,r,mu1 :
 option remember:

  if (m<0 or n<0 or m<n) then RETURN({}): fi:

  if m=0 then RETURN({}): fi:

  gu:={}:

  for r from 1 to m do 
    mu:=SetU(m-r,n) union SetD(m-r,n): 
    gu:= gu union {seq([op(mu1),0$r],mu1 in mu)}: 
  od:

  gu:

  end:


#u(m,n): inputs non-negative integers m,n
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step.
u:=proc(m,n) local r:
  option remember:
  if m=0 then 0:
  else add(d(m-r,n-r)+f(m-r,n-r),r=1..n):
fi:
end:


#d(m,n): inputs non-negative integers m,n
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step.
d:=proc(m,n) local r: 
option remember:
  if m=0 then 
    if n=0 then RETURN(1):
     else RETURN(0):
     fi:
   fi:
 add(u(m-r,n+r)+f(m-r,n+r),r=1..m):
end:

#f(m,n): inputs non-negative integers m,n
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step.
f:=proc(m,n) local r:
option remember:
  if m=0 then RETURN(0):
   fi:
    add(u(m-r,n)+d(m-r,n),r=1..m):
end:


#Ru(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step
#s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

Ru:=proc(m,n,A,B,F) local gu,mu,r,mu1 :
  option remember:
  if (m<0 or n<0 or m<n) then RETURN({}): fi:
  if m=0 and n=0 then RETURN({}): fi:
  gu:={}:

  if F={} or F=d then
    if member(n,A) then RETURN({}): fi: fi:

  for r from 1 to n do 
    mu:=Rd(m-r,n-r,A,B,u) union Rf(m-r,n-r,A,B,u):
    gu:= gu union {seq([op(mu1),1$r],mu1 in mu)}:  
    od:
  gu:
end:

#Rd(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers, and F={} or d or f
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step
#s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

Rd:=proc(m,n,A,B,F) local gu,mu,r,mu1 :
 option remember:
  if (m<0 or n<0 or m<n) then RETURN({}): fi:
  if m=0 then 
    if n=0 then RETURN({[]}):
    else RETURN({}): fi: fi:

  gu:={}:

  if F={} or F=u then
    if member(n,B) then RETURN({}): fi: fi:

  for r from 1 to m do 
    mu:=Ru(m-r,n+r,A,B,d) union Rf(m-r,n+r,A,B,d):
    gu:= gu union {seq([op(mu1),-1$r],mu1 in mu)}: 
  od:
  
  gu:
end:


#Rf(m,n,A,B,F): inputs non-negative integers m and n, finite sets A and B of non-negative integers , and F={} or d or f
#Outputs the set of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step
#s.t. no peak heights belongs to A and no valley height belongs to B when F is appended to the walk.
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

Rf:=proc(m,n,A,B,F) local gu,mu,r,mu1, nu, gu1, gu2 :
 option remember:
  if (m<0 or n<0 or m<n) then RETURN({}): fi:
  if m=0 then RETURN({}): fi:

  if member(n,B) and n=0 then RETURN({[0$m]}):
  elif F={} and member(n,B) and member(n,A) then RETURN({}):
  else
    gu1:={}: gu2:={}: 

    for r from 1 to m do 
      mu:=Ru(m-r,n,A,B,f):  nu:=Rd(m-r,n,A,B,f):
      gu1:= gu1 union {seq([op(mu1),0$r],mu1 in mu)}:
      gu2:= gu2 union {seq([op(nu1),0$r],nu1 in nu)}:
    od:

    if F=u and member(n,B) then gu:=gu1:
    elif F=d and member(n,A) then gu:=gu2:
    elif F={} and member(n,B) then gu:=gu1:
    elif F={} and member(n,A) then gu:=gu2:
    else
     gu:=gu1 union gu2:
    fi:
 fi:
gu:
end:


#NRu(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or d or f
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step
#s.t. no peak has height belongs to A and no valley height belongs to B
#no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRu:=proc(m,n,A,B,C,D,E,F) local r, gu:
  option remember:

  if (m<0 or n<0 or m<n) then RETURN(0): fi:
  if m=0 then RETURN(0): fi:

  gu:=0:
  if F={} or F=d then
    if member(n,A) then  RETURN(0): fi: fi:

  for r from 1 to n do 
    if not member(r,C) then 
    	gu:=gu+NRd(m-r,n-r,A,B,C,D,E,u)+NRf(m-r,n-r,A,B,C,D,E,u): fi: od:
  gu:

  end:

#NRd(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers, and F={} or u or f
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step
#s.t. no peak has height belongs to A and no valley height belongs to B
#no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRd:=proc(m,n,A,B,C,D,E,F) local r, gu:
  option remember:

  if (m<0 or n<0 or m<n) then RETURN(0): fi:
  if m=0 then 
    if n=0 then RETURN(1):
    else RETURN(0): fi: fi:
   gu:=0:

  if F={} or F=u then
    if member(n,B) then RETURN(0): fi: fi:

  for r from 1 to m do 
  	if not member(r,D) then  
       gu:= gu + NRu(m-r,n+r,A,B,C,D,E,d)+ NRf(m-r,n+r,A,B,C,D,E,d): 
     fi: od:
  gu:
  end:

#NRf(m,n,A,B,C,D,E,F): inputs non-negative integers m and n, finite sets A, B,C,D, and E of non-negative integers , and F={} or u or d
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step
#s.t. no peak has height belongs to A and no valley height belongs to B
#no upward run belongs to C, no downward run belongs to D, and no flat run belongs to E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRf:=proc(m,n,A,B,C,D,E,F) local r, gu,gu1,gu2:
 option remember:

  if (m<0 or n<0 or m<n) then RETURN(0): fi:

  if m=0 then RETURN(0): fi:

  if member(n,B) and n=0 then gu:=1:
  elif F={} and member(n,B) and member(n,A) then gu:=0:
  else
    gu1:=0: gu2:=0: 

    for r from 1 to m do 
      if not member(r,E) then  
        gu1:= gu1 + NRu(m-r,n,A,B,C,D,E,f): 
        gu2:=  gu2 + NRd(m-r,n,A,B,C,D,E,f): 
    fi:od:

    if F=u and member(n,B) then gu:=gu1:
    elif F=d and member(n,A) then gu:=gu2:
    elif F={} and member(n,B) then gu:=gu1:
    elif F={} and member(n,A) then gu:=gu2:
    else gu:=gu1+gu2:
  fi: fi:

  gu:
end:


#IsRange1(f,k,n): inputs a linear expression f in k and a non-negative integer n
#decides whether f(k0)=n for some non-negative integer k0.
#This was taken from Dr. Doron Zeilberger's Dyck.txt

IsRange1:=proc(f,k,n) local k0:
  k0:=solve(f=n,k):

  if k0>=0 and type(k0,integer) then
    true:
  else
   false:
  fi:
end:

#IsRange(F,k,n): inputs a set of linear expressions F={f} in k and a non-negative integer n
#decides whether f(k0)=n for some non-negative integer k0 and some f in F
#This was taken from Dr. Doron Zeilberger's Dyck.txt

IsRange:=proc(F,k,n) local f:
  for f in F do
    if IsRange1(f,k,n) then
      RETURN(true):
    fi:
  od:
  false:
end:


#NRur(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or d or f
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in an up step
#s.t. no peak has height is in the range of A and no valley height is in the range of B
#no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRur:=proc(m,n,A,B,C,D,E,F,r) local gu,k:
  option remember:

  if (m<0 or n<0 or m<n) then
    RETURN(0):
  fi:

  if m=0 and n=0 then
    RETURN(0):
  fi:

  if F={} or F=d then 
    if IsRange(A,r,n) then
      RETURN(0):
  fi: fi:

  gu:=0:

  for k from 1 to n do
    if not IsRange(C,r,k) then
      gu:=gu+NRdr(m-k,n-k,A,B,C,D,E,u,r)+NRfr(m-k,n-k,A,B,C,D,E,u,r):
    fi:
  od:

  gu:
end:


#NRdr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or f
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a down step
#s.t. no peak has height is in the range of A and no valley height is in the range of B
#no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRdr:=proc(m,n,A,B,C,D,E,F,r) local gu,k:
  option remember:

  if (m<0 or n<0 or m<n) then
    RETURN(0):
  fi:

  if m=0 and n=0 then
    RETURN(1):
  fi:

  if F={} or F=u then
    if IsRange(B,r,n) then
      RETURN(0):
  fi: fi:

  gu:=0:

  for k from 1 to m do
    if not IsRange(D,r,k) then
      gu:=gu+NRur(m-k,n+k,A,B,C,D,E,d,r)+NRfr(m-k,n+k,A,B,C,D,E,d,r):
   fi:
  od:

  gu:
end:

#NRfr(m,n,A,B,C,D,E,F,r): inputs non-negative integers m and n, sets A,B,C,D,E of linear expressions in the variable r, and F={} or u or d
#Outputs the number of walks from (0,0) to (m,n) which stay weakly above the x-axis and end in a flat step
#s.t. no peak has height is in the range of A and no valley height is in the range of B
#no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

NRfr:=proc(m,n,A,B,C,D,E,F,r) local gu,k,gu1,gu2:
  option remember:

  if (m<0 or n<0 or m<n) then
    RETURN(0):
  fi:

  if m=0 and n=0 then
    RETURN(0):
  fi:

  if IsRange(B,r,n) and n=0 then gu:=1:
  elif F={} and IsRange(B,r,n) and IsRange(A,r,n) then gu:=0:
  else
    gu1:=0: gu2:=0: 

    for k from 1 to m do
      if not IsRange(E,r,k) then
        gu1:= gu1 + NRur(m-k,n,A,B,C,D,E,f,r):  
        gu2:=  gu2 + NRdr(m-k,n,A,B,C,D,E,f,r):
    fi:od:
 
    if F=u and IsRange(B,r,n) then gu:=gu1: 
    elif F=d and IsRange(A,r,n) then gu:=gu2:
    elif F={} and IsRange(B,r,n) then gu:=gu1:
    elif F={} and IsRange(A,r,n) then gu:=gu2:
    else gu:=gu1+gu2:
    fi:
  fi:
  gu:
end:


#SeqABCDE(A,B,C,D,E,K): inputs finite sets A,B,C,D,E of non-negative integers and positive integer K
#Outputs the first K terms of the sequence of Motzkin paths of length n
#where no peak height belongs to A and no valley height belongs to B
#No upward run belongs to C, No downward run in D, and no flat run belongs to E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk
SeqABCDE:=proc(A,B,C,D,E,K) local i,L,m:
  L:=[seq(NRd(m,0,A,B,C,D,E,{})+NRf(m,0,A,B,C,D,E,{}),m=0..K)]:
end:



#SeqABCDEr(A,B,C,D,E,r,K): inputs sets A,B,C,D,E of linear expressions in the variable r and positive integer K
#Outputs the first K terms of the sequence of Motzkin paths of length n
#where no peak has height is in the range of A and no valley height is in the range of B
#no upward run is in the range of C, no downward run is in the range of D, and no flat run is in the range of E
#Here, a peak is defined as any chain U F^k D where F^k is a flat run of length k where k is a non-negative integer
#Similarly, a valley is a chain D F^k U where F^k is a flat run of length k where k is a non-negative integer
#or a chain D F^k that occurs at the end of a walk

SeqABCDEr:=proc(A,B,C,D,E,r,K) local i,L:

L:=[seq(NRdr(m,0,A,B,C,D,E,{},r)+NRfr(m,0,A,B,C,D,E,{},r),m=0..K)]:

end:



#####Start Brute force#####
 
#Motzkin(n): outputs all Motzkin paths of length n.
Motzkin:=proc(n) local k:
  option remember:
  if n<0 then {}:
  elif n=0 then {[]}:
  else
    {seq(op(Mk(n,k)),k=1..n)}: 
  fi:
end:



#Mk(n,k): All Motzkin paths of length n of the form 1 L1 (-1) L3 or L1 0 L2  
#where L1 in Motzkin(k-1), L2 in Motzkin(n-k), and L3 ni Motzkin (k-2).

Mk:=proc(n,k) local gu1,gu2,mu1,mu2,gu3, mu3:
option remember:
  if n<0 or n=0 then RETURN({}): fi:

  if k=0 then RETURN({}): fi:
  
  gu1:=Motzkin(k-1):
  gu2:=Motzkin(n-k):
  gu3:=Motzkin(k-2):

  {seq(seq([1, op(mu3), -1, op(mu2)], mu3 in gu3), mu2 in gu2),seq(seq([op(mu1), 0, op(mu2)], mu1 in gu1), mu2 in gu2)}:
end:



#Targem1(L): translates a Motzkin path into a sequence of positive, negative integers, and sets
#representing the upwards, downwards, and flat runs respectively
#such that the partial sums of the integers are always non-negative

Targem1:=proc(L) local i,j,k1,k2,k3:
  if L=[] then RETURN([]): fi:

  for k1 from 1 to nops(L) while L[k1]=0 do od:

  for i from k1 to nops(L) while L[i]=1 do od:

  for k2 from i to nops(L) while L[k2]=0 do od:
  
  for j from k2 to nops(L) while L[j]=-1  do od:

  for k3 from j to nops(L) while L[k3]=0 do od:
  subs({0}=NULL,0=NULL,[{k1-1},i-k1,{k2-i},k2-j, {k3-j}, op(Targem1([op(k3..nops(L),L)]))]):

end:



#IsGoodCDE(L,C,D,E): Inputs a Motzkin path L, and finite sets C,D,E of non-negative integers.
#Decides whether no upwards run in L belongs to C, no downwards run in L belongs to D, and no flat run in L belongs to E.
IsGoodCDE:=proc(L,C,D,E) local L1,i:
  L1:=Targem1(L):

  for i from 1 to nops(L1) do
    if type(L1[i],set)=true then 
      if member(op(L1[i]),E) then RETURN(false): fi:
    elif type(L1[i],positive)=true then 
      if member(L1[i],C) then RETURN(false): fi:
    else
      if member(-L1[i],D) then RETURN(false): fi:
    fi:
  od:

  true:
end:


#IsGoodCDEr(L,C,D,E,r): Inputs a Motzkin path L, and sets C,D,E of linear expressions in the variable r.
#Decides whether no upwards run in L belongs to the range of C, no downward run in L belongs to the range of D,
#and no flat run in L belongs to the range of E.

IsGoodCDEr:=proc(L,C,D,E,r) local L1,i:
  L1:=Targem1(L):

 for i from 1 to nops(L1) do
   if type(L1[i],set)=true then 
     if IsRange(E,r,op(L1[i])) then RETURN(false): fi:
   elif type(L1[i],positive)=true then 
     if IsRange(C,r,L1[i]) then RETURN(false): fi:
   else
    if IsRange(D,r,-L1[i]) then RETURN(false): fi:
   fi:
  od:

  true:
end:


#Targem2(L): translates a Motzkin path into a sequence of the partial sums when it between a 
#chain with steps in {0,1} starting with 1 (i.e. weakly increasing chain starting with an upwards step)
#and a chain with steps in {0,-1} starting with -1 (i.e. a weakly decreasing chain starting with a downwards step)

Targem2:=proc(L) local L1,i,L2,su:
  L1:=Targem1(L):
  L2:=[]:
  su:=0:
  
  for i from 1 to nops(L1) do 
    if type(L1[i],set)=false then 
      su:=su+L1[i]: 
      L2:=[op(L2),su]:
    else
      if i < nops(L1) and i>1 then
        if type(L1[i-1],positive)=type(L1[i+1],positive) then
        	su:=su+L1[i+1]:
        	L2:=[op(1..-2,L2),su]: 
        	i:=i+1: 
  fi: fi: fi: od:
  L2:
end:


#IsGoodA(L,A): Inputs a Motzkin path L, and a finite set A of non-negative integers.
#Decides whether none of the heights in a peak location of L belongs to A
#A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1,
#s.t. any steps between (if any) are 0.

IsGoodA:=proc(L,A) local L1,i:
  L1:=Targem2(L):

  for i from 1 to nops(L1)/2 do
    if member(L1[2*i-1],A) then
      RETURN(false):
    fi:
  od:
  
  true:
end:


#IsGoodB(L,B): Inputs a Motzkin path L, and a finite set B of non-negative integers,
#Decides whether none of the heights in a valley location of L belongs to B.
#A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1,
#s.t. any steps between (if any) are 0.
IsGoodB:=proc(L,B) local L1,i:
  L1:=Targem2(L):

  for i from 1 to nops(L1)/2 do
  	if member(L1[2*i],B) then RETURN(false): fi:
  od:

  true:
end:



#IsGoodAr(L,A,r): Inputs a Motzkin path L, and a set A of linear expressions in the variable r.
#Decides whether none of the heights in a peak location belongs to the range of A
#A peak height is the maximal height of a chain beginning with the step 1 and ending with the step -1,
#s.t. any steps between (if any) are 0.

IsGoodAr:=proc(L,A,r) local L1,i:
  L1:=Targem2(L):

  for i from 1 to nops(L1)/2 do
    if IsRange(A,r,L1[2*i-1]) then
      RETURN(false):
    fi:
  od:

  true:
end:


#IsGoodBr(L,A,r): Inputs a Motzkin path L, and a set B of linear expressions in the variable r. 
#Decides whether none of the heights in a valley location belongs to the range of B
#A valley height is the minimal height of a chain beginning with the step -1 and ending with the step 1,
#s.t. any steps between (if any) are 0.

IsGoodBr:=proc(L,B,r) local L1,i:
  L1:=Targem2(L):

  for i from 1 to nops(L1)/2 do
    if IsRange(B,r,L1[2*i]) then
      RETURN(false):
    fi:
  od:

  true:
end:

#MotzkinABCDE(A,B,C,D,E,n): Inputs five finite sets A, B, C,D,E of non-negative integers and a positive integer n,
#Outputs the set of Motzkin paths counted by SeqABCDE(A,B,C,D,E,K). In other words:
#where no peak can be of height that belongs to A and no valley can be of height that belongs to B
#No length of upward run belongs to C and no length of downward run belongs to D and no length of flat run belongs to E
#A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down),
#s.t. any steps between (if any) are 0 (i.e. flat).

MotzkinABCDE:=proc(A,B,C,D,E,n) local gu,mu,mu1:
  mu:=Motzkin(n):
  gu:={}:
  for mu1 in mu do
    if (IsGoodCDE(mu1,C,D,E) and IsGoodA(mu1,A) and IsGoodB(mu1,B)) then
      gu:=gu union {mu1}:
    fi:
  od:

  gu:
end:



#MotzkinABCDEr(A,B,C,D,E,r,n): Inputs five sets  A, B, C,D,E of linear expressions in the variable r, and a positive integer n.
#Outputs the set of Motzkin paths counted by SeqABCDEr(A,B,C,D,E,r, K). In other words:
#where no peak can be of height that belongs to the range of A and no valley can be of height that belongs to  the range of B,
#No length of upward run belongs to the range of C and no length of downward run belongs to the range of D
#and no length of flat run belongs to the range of E.
#A peak height is the maximal height of a chain beginning with the step 1 (i.e. up) and ending with the step -1 (i.e. down),
#s.t. any steps between (if any) are 0 (i.e. flat).

MotzkinABCDEr:=proc(A,B,C,D,E,r,n) local gu,mu,mu1:
  mu:=Motzkin(n):
  gu:={}:
  for mu1 in mu do
    if (IsGoodCDEr(mu1,C,D,E,r)  and IsGoodAr(mu1,A,r) and IsGoodBr(mu1,B,r)) then
      gu:=gu union {mu1}:
    fi:
  od:

  gu:
end:


#SeqABCDEbrute(A,B,C,D,E,K): Same output as SeqABCDE(A,B,C,D,E,K) but by complete brute force. 
#FOR CHECKING PURPOSES ONLY
#Don't  make K too big. 
SeqABCDEbrute:=proc(A,B,C,D,E,K) local i:
  [seq(nops(MotzkinABCDE(A,B,C,D,E,i)),i=0..K)]:
end:

#SeqABCDErBrute(A,B,C,D,K):  Same output as SeqABCDEr(A,B,C,D,E,r,K) but by complete brute force. 
#FOR CHECKING PURPOSES ONLY
#Don't  make K too big. 

SeqABCDErBrute:=proc(A,B,C,D,E,r,K) local i:
[seq(nops(MotzkinABCDEr(A,B,C,D,E,r,i)),i=0..K)]:

end:





###Start from SCHUTZENBERGER.txt
#empir(gu,degx,degP,x,P) 
#to "fit" an algebraic equation F(P(x),x)=0 of degree
#degP in P(x) and degx in n for P(x):=sum_i gu[i]*x^i
 
empir:=proc(gu,degx,degP,x,P)
local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip,ka,mu1,Ftry:

if (1+degx)*(1+degP) > nops(gu)-15 then
 RETURN(`sequence too small`):
fi:
 
F:=0:
var:={}:
for i1 from 0 to degx do
 for i2 from 0 to degP do
  F:=F+a[i1,i2]*x^i1*P^i2:
  var:=var union {a[i1,i2]}:
 od:
od:
 
cand:=0:
 
for i1 from 0 to nops(gu)-1 do
 cand:=cand+op(i1+1,gu)*x^i1:
od:
 
lu:=subs(P=cand,F):
lu:=taylor(lu,x=0,nops(gu)-1):
 
eq:={}:
 
for i1 from 0 to nops(gu)-2 do
eq:=eq union {coeff(lu,x,i1)=0}
od:
 
mu:=solve(eq,var):

ka:=0:

for mu1 in mu do
 if op(1,mu1)=op(2,mu1) then
 ka:=ka+1:
 fi:
od:

if ka>1 then
 RETURN(FAIL):
fi:

F:=subs(mu,F):
 
if F=0 then
 RETURN(FAIL):
fi:
 
flo:=degree(F,P):
pip:=coeff(F,P,flo):
flo:=degree(pip,x):
pip:=coeff(pip,x,flo):
F:=normal(F/pip):

Ftry:=taylor(subs(P=add(gu[i1+1]*x^i1,i1=0..nops(gu)-1),F),x=0,nops(gu)+1):

if {seq(coeff(Ftry,x,i1),i1=0..nops(gu)-2)}<>{0} then
 RETURN(FAIL):
fi:

normal(F):
end:



Aeq:=proc(gu,x,P)
local degx,degP,L,lu:
 
for L from 1 to (nops(gu)-15)/3 do
for degP from 1 to L do
for degx from 0 to min(trunc((nops(gu)-15)/(1+degP))-1,L-degP) do
 
lu:=empir(gu,degx,degP,x,P):
 
if lu<>FAIL then
RETURN(lu):
fi:
od:
od:
od:
FAIL:
end:
 
###end from SCHUTZENBERGER.txt
 
   
 
##start from Findrec
ezraFindrec:=proc()
if args=NULL then

print(` FindRec.txt: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of ONE Variable`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`  findrec, Findrec, FindrecF`):
print():

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):


fi:
 
end:


###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
#ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
RETURN(FAIL):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
#ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
RETURN(FAIL):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
ope:=ope[1]:

if {seq(add(subs(n=n1+i1,coeff(ope,N,i1))*f[n1+i1],i1=0..degree(ope,N)),n1=1..nops(f)-degree(ope,N))}<>{0} then
 RETURN(FAIL):
fi:

 ope:=Yafe(ope,N)[2]:


if SeqFromRec(ope,n,N,[op(1..degree(ope,N),f)],nops(f))<>f then
 RETURN(FAIL):
else
RETURN(ope):
fi:
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
# print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL and degree(ope,N)<>0 then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1,ka,i1:
ope1:=ope:

L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ka:={seq(subs(n=i1,lcoeff(ope1,N)),i1=1..K)}:

if member(0,ka) then
 RETURN(FAIL):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
if member(0, { seq( subs(n=n1,denom(coeff(ope1,N,-j1)  ))  ,j1=1..L )}) then
 RETURN(FAIL):
else
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
fi:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
#ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
RETURN(FAIL):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:
###End from Findrec